Introduction to Bayesian pharmacometric data analysis using NONMEM
PopPK analysis using MCMC
Metrum Research Group
October 15, 2019
Population Parameter Estimation
Estimation of pediatric atorvastatin PK given informative priors based on adults
Specification of prior distributions for population means and inter-individual covariance
Let’s revise the priors to be more consistent with current practice. Specifically,
- Let’s use weakly informative priors to guard against implausible values and help improve sampling performance.
- Update variance of CL/F, V2/F, and Q/F to 2.
- Add variance of \(F_{chew}\) and SD to 1.
- Let’s use a LKJ prior for IIV correlation matrix and lognormal priors for diag(OMEGA).
- OLKJDF = 1 and OVARF=1.
Setup R environment
rm(list = ls())
gc()
modelName <- "atorv1"
scriptName <- paste(modelName, "Rmd", sep = ".")
fitModel <- FALSE
## Relative paths assuming the working directory is the script directory
## containing this script
scriptDir <- getwd()
projectDir <- dirname(scriptDir)
dataDir <- file.path(projectDir, "data")
modelDir <- file.path(projectDir, "model")
outDir <- file.path(modelDir, modelName)
figDir <- file.path(projectDir, "deliv", "figure", modelName)
tabDir <- file.path(projectDir, "deliv", "table", modelName)
invisible(dir.create(figDir, recursive = TRUE))## Warning in dir.create(figDir, recursive = TRUE): '/data/IntroBayesNONMEM/
## deliv/figure/atorv1' already existsinvisible(dir.create(tabDir, recursive = TRUE))## Warning in dir.create(tabDir, recursive = TRUE): '/data/IntroBayesNONMEM/
## deliv/table/atorv1' already exists.libPaths("lib")
library(metrumrg)## Loading required package: reshape## Loading required package: lattice## Loading required package: XML## Loading required package: MASS## Loading required package: grid## metrumrg 5.57## enter "?metrumrg" for helpsuppressMessages(library(rstan))
suppressMessages(library(bayesplot))
suppressMessages(library(tidyverse))
library(gridExtra)##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combinelibrary(parallel)
library(knitr)
library(kableExtra)##
## Attaching package: 'kableExtra'## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(qapply)## Loading required package: rlecuyer##
## Attaching package: 'qapply'## The following object is masked from 'package:metrumrg':
##
## qstatlibrary(PKPDmisc)
knitr::opts_chunk$set(
comment = '.',
fig.height = 5,
fig.width = 9,
eval.after = 'fig.cap',
message = FALSE,
fig.path = file.path(figDir, paste(modelName, "_", sep = ""))
)
## Go back to default ggplot2 theme that was overridden by bayesplot
theme_set(theme_gray())
source(file.path(scriptDir, "tools", "stanTools2.R"))
source(file.path(scriptDir, "tools", "functions.R"))
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
set.seed(11191951)
set.seed(10271998) The NONMEM model stub
modelFile <- file.path(modelDir, paste(modelName, "stub.ctl", sep = ""))
modelText <- readLines(modelFile)
cat(modelText, sep = "\n"). $PROBLEM RUN# atorv1
. $INPUT C ID TIME AMT DV CMT DOSE EVID MDV ADDL II SS FORM WT PER BLQ
. ;$INPUT C ID SID=DROP STDY=DROP TIME DV DVR=DROP DVC=DROP AMT AMTO=DROP CMT TYPE=DROP
. ; DOSE DOSN=DROP NTIM=DROP EVID MDV ADDL II SS FORM WT SEX=DROP AGE
. ; HT=DROP RACE=DROP PER BLQ TRCD=DROP ANLT=DROP TXT=DROP
. ;$DATA ../../../data/atorv1172atvb.csv IGNORE=(C='C', BLQ.EQ.1)
. $DATA ../../../data/atorvWrkShop.csv IGNORE=(C='C', BLQ.EQ.1)
.
. $ABBR DECLARE INTEGER FIRST_WRITE_PAR
. $ABBR DECLARE INTEGER FIRST_WRITE_IPAR
.
. $SUBROUTINES ADVAN4 TRANS4 OTHER = ../../extrasend.f90
.
. $PK
. ; Request extra information for Bayesian analysis.
. ; An extra call will then be made for accepted samples
. include '/opt/NONMEM/nm74/util/nonmem_reserved_general'
. BAYES_EXTRA_REQUEST=1
.
. nThin = THETA(8)
.
. VWT = LOG(WT/70) ; normalized to 70 kg adult
. MU_1 = THETA(1) + VWT * 0.75 ; CL
. MU_2 = THETA(2) + VWT ; V2
. MU_3 = THETA(3) + VWT * 0.75 ; Q
. MU_4 = THETA(4) + VWT ; V3
. MU_5 = THETA(5) ; ka
. MU_6 = THETA(6) ; F1chew
. MU_7 = THETA(7) ; SD
.
. " CALL EXTRASEND()
.
. ;ATORVASTATIN
.
. CL = EXP(MU_1 + ETA(1)) ; atv cl
. V2 = EXP(MU_2 + ETA(2))
. Q = EXP(MU_3 + ETA(3))
. V3 = EXP(MU_4 + ETA(4))
. ;deltaKA = EXP(MU_5 + ETA(5))
. KA = EXP(MU_5 + ETA(5))
. F1chew = EXP(MU_6 + ETA(6))
. SD = EXP(MU_7 + ETA(7))
.
. ; Constrain KA > lambda1 (typical value) to enhance identifiability
. ;T1 = CL/V2
. ;T23 = Q/V2
. ;T32 = Q/V3
. ;L1 = ((T1+T23+T32)+SQRT((T1+T23+T32)**2-4*T1*T32))/2
. ;L2 = ((T1+T23+T32)-SQRT((T1+T23+T32)**2-4*T1*T32))/2
.
. ;KA = deltaKA + L2
. S2 = V2/1000 ; DOSE IN uM & CONC IN nM/L
. F1 = 1.0 ; TABLET AS REFERENCE
. IF(FORM.EQ.2) F1 = F1chew ; CHEWABLE F1
.
. IF(BAYES_EXTRA==1 .AND. NIREC==1 .AND. NDREC==1 .AND. &
. ITER_REPORT>0 .AND. &
. ITER_REPORT/nThin == INT(ITER_REPORT/nThin)) THEN
. IF(FIRST_WRITE_PAR==0) THEN
. " OPEN(unit=52,FILE='./par.txt')
. FIRST_WRITE_PAR=1
. ENDIF
. " WRITE(52,'(I12,1X,20(1X,1PG19.10E3))') ITER_REPORT, &
. " OMEGA(1,1), OMEGA(2,1), OMEGA(2,2), &
. " THETA(1), THETA(2), THETA(3), THETA(4), THETA(5), THETA(6), THETA(7)
. ENDIF
. IF(BAYES_EXTRA==1 .AND. NDREC==1 .AND. ITER_REPORT>0 .AND. &
. ITER_REPORT/nThin == INT(ITER_REPORT/nThin)) THEN
. IF(FIRST_WRITE_IPAR==0) THEN
. " OPEN(unit=50,FILE='./ipar'//TFI(PNM_NODE_NUMBER)//'.txt')
. FIRST_WRITE_IPAR=1
. ENDIF
. " WRITE(50,'(I12,1X,F14.0,10(1X,1PG19.10E3))') ITER_REPORT, ID, &
. " ETA(1), ETA(2)
. ENDIF
.
.
. $ERROR
.
. LOQ = 0.45 ; nM/L
. IPRED = F
. DUM = (LOQ-IPRED)/(SD*IPRED)
. CUMD=PHI(DUM)
. IF(BLQ.EQ.0)THEN
. F_FLAG=0
. Y = F*(1+SD*ERR(1)) ; ATV
. ELSE
. F_FLAG=1
. Y=CUMD
. ENDIF
.
. ; Code below into chainsAssemble initial estimates, priors and table specs
## create initial estimates
geninit <- function(){
paste(c(
"; Initial values of THETA",
"$THETA",
paste(rnorm(1, log(1000), 0.5), " ; TOTAL CL ATV = THETA(1)"),
paste(rnorm(1, log(600), 0.5), " ; V2 = THETA(2)"),
paste(rnorm(1, log(130), 0.5), " ; Q = THETA(3)"),
paste(rnorm(1, log(1500), 0.5), " ; V3 = THETA(4) "),
paste(rnorm(1, log(1.2), 0.5), " ; KA= THETA(5) "),
paste(rnorm(1, log(1), 0.25), " ; THETA(6) relative F1 CHEWABLE = THETA(5) "),
paste(log(0.2), " ; THETA(7) - PRO RES ERROR"),
"$OMEGA BLOCK(2) ;INITIAL values of OMEGAs",
"0.1003 ; cl",
".06034 0.8477 ; v2",
"$OMEGA 0 FIX",
"$OMEGA 0 FIX",
"$OMEGA 0 FIX",
"$OMEGA 0 FIX",
"$OMEGA 0 FIX",
"$SIGMA ;Initial value of SIGMA",
"1 FIX; pro"),
sep = "\n")
}
## Set parameters of the prior distribution
priors <- paste(c("$PRIOR NWPRI",
"$THETAP ; Prior information of THETAS",
paste("(", 5.65, " FIX) ; THETA(1) TOTAL CL ATV"),
paste("(", 5.60, " FIX) ; THETA(2) V2"),
paste("(", 4.29, " FIX) ; THETA(3) Q"),
paste("(", 7.07, " FIX) ; THETA(4) V3"),
paste("(", -1.33, " FIX) ; THETA(5) KA"),
paste("(", log(1), " FIX) ; THETA(6) F1 chewable"),
paste("(", log(0.2), " FIX) ; THETA(7) log(sigma)"),
"$THETAPV BLOCK(7) ; variances for priors on THETAS (var-cov)",
"2 FIX ; CL weakly informative",
"0.00 2 ; V2 weakly informative",
"0.00 0.00 2 ; Q weakly informative",
"0.00 0.00 0.00 .0156 ; V3 SE**2 informative",
"0.00 0.00 0.00 0.00 .00684 ; KA SE**2 informative",
"0.00 0.00 0.00 0.00 0.00 1 ; F1 chewable weakly informative",
"0.00 0.00 0.00 0.00 0.00 0.00 1 ; log(sigma) weakly informative",
## Parameters LKJ / lognormal prior; off-diagonal and OMEGAP ignored
"$OMEGAP BLOCK(2) ; prior information for OMEGA",
".1003 FIXED ; cl",
".06034 .8477 ; v2" ,
"$OMEGAPD (4.0, FIXED) ; df for OMEGA prior"
## Parameters for inverse Wishart prior
## "$OMEGAP BLOCK(2) ; prior information for OMEGA",
## ".1003 FIXED ; cl",
## ".06034 .8477 ; v2" ,
## "$OMEGAPD (4.0, FIXED) ; df for OMEGA prior"
),
sep = "\n")
## Table statements
tables <- paste(c("$TABLE ID EVID TIME DV IPRED CWRES CWRESI NPDE WT",
"NOPRINT ONEHEADER FILE=./.tab",
"$TABLE ID WT CL V2 Q KA V3 ETA1 ETA2 PER FORM",
"NOPRINT ONEHEADER FILE=./par.tab"),
sep = "\n")
##tables <- paste(c("$TABLE ID EVID TIME DV IPRED CWRES CWRESI NPDE AGE WT",
## "NOPRINT ONEHEADER FILE=./.tab",
## "$TABLE ID WT CL V2 Q KA V3 ETA1 ETA2 PER FORM",
## "NOPRINT ONEHEADER FILE=./par.tab"),
## sep = "\n")Run NONMEM
nChains <- 4
nPost <- 250 ## Number of post-burn-in samples per chain after thinning
nBurn <- 250 ## Number of burn-in samples per chain after thinning
nThin <- 1
seed = sample(10000, size = nChains)
if(fitModel){
if(!file.exists(file.path(modelDir, modelName)))
dir.create(file.path(modelDir, modelName))
runChains <- mclapply(1:nChains, runChain,
modelName = modelName, modelDir = modelDir,
priors = priors, tables = tables,
nPost = nPost, nBurn = nBurn, nThin = nThin,
seed = seed, print = 100,
## LKJ prior for IIV correlation matrix with oarameter = OLKJDF
OLKJDF = 1,
## lognormal prior for IIV SDs
## log(sqrt(Omega(i))) ~ Normal(log(sqrt(OmegaPrior(i))), 1/OVARF)
OVARF = 1, AUTO = 2,
NUTS_DELTA = 0.8,
grid = FALSE,
method = "NUTS",
## pe = "orte 16",
mode = "nonpara")
# NONR(run = paste(modelName, ".", 1, sep = ""),
# command = "/opt/NONMEM/nm74/nmqual/autolog.pl",
# project = file.path(modelDir, modelName),
# grid = FALSE,
# wait = FALSE,
# diag = FALSE,
# fdata = FALSE,
# purge = FALSE,
# checkrunno = TRUE)
}read in data file
dataFile <- "atorvWrkShop.csv"
nmData <- read.csv(file.path(dataDir, dataFile), as.is = TRUE)Get population parameters
chains <- 1:nChains
##chains <- c(1, 4)
nChains2 <- length(chains)
## Read in posterior distributions by chain and add a column called chain
popParameters <- map(chains, function(thisChain){
param <- data.table::fread(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = ""),
"par.txt"), data.table = FALSE)
names(param) <- c("iteration", "OMEGA11", "OMEGA21", "OMEGA22",
"THETA1", "THETA2", "THETA3", "THETA4", "THETA5",
"THETA6", "THETA7")
param %>% mutate(chain = thisChain)
}) %>% bind_rows %>%
mutate(sample = 1:n())
## Calculate more interpretable paraneters
popParameters2 <- popParameters %>%
mutate(CLHat = exp(THETA1),
V2Hat = exp(THETA2),
QHat = exp(THETA3),
V3Hat = exp(THETA4),
## k10 = CLHat / V2Hat,
## k12 = QHat / V2Hat,
## k21 = QHat / V3Hat,
## lambda = ((k10 + k12 + k21) -
## sqrt((k10 + k12 + k21)^2 - 4 * k10 * k21)) / 2,
## KAHat = exp(THETA4) + lambda,
## KAHat = exp(THETA4 + lambda),
KAHat = exp(THETA5),
F1Chewable = exp(THETA6),
sigma1 = exp(THETA7)) %>%
select(starts_with("OMEGA"),
CLHat, V2Hat, QHat, V3Hat, KAHat, F1Chewable, sigma1, chain, iteration, sample)
## Convert to 3-D array with dims = {iterations, chains, parameters}.
popParArray <- array(double(nPost * nChains2 * (ncol(popParameters2) - 3)),
dim = c(nPost, nChains2, ncol(popParameters2) - 3),
dimnames = list(NULL, NULL, setdiff(names(popParameters2), c("chain", "iteration", "sample"))))
for(iChain in 1:nChains2){
popParArray[,iChain,] <- popParameters2 %>%
filter(chain == chains[iChain],
iteration > 0) %>%
select(-chain, -iteration, -sample) %>%
as.matrix
}Get individual parameters
indParameters = map(chains, function(thisChain){
files <- list.files(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = "")))
iparFiles <- files[grepl("^ipar", files)]
param <- map(iparFiles, function(thisFile){
data.table::fread(file.path(modelDir, modelName, paste(modelName, ".", thisChain, sep = ""),
thisFile), data.table = FALSE)
}) %>%
bind_rows
names(param) <- c("iteration", "ID", "ETA1", "ETA2")
param %>% mutate(chain = thisChain)
}) %>% bind_rows %>%
arrange(chain, iteration, ID) %>%
left_join(popParameters)MCMC diagnostics and posterior distributions of parameters
options(bayesplot.base_size = 12,
bayesplot.base_family = "sans")
color_scheme_set(scheme = "brightblue")
myTheme <- theme(text = element_text(size = 12),
axis.text = element_text(size = 12))
plot_mcmcHistory <- mcmcHistory(popParArray,
nParPerPage = 5, myTheme = myTheme)
plot_mcmcDensityByChain <- mcmcDensity(popParArray,
nParPerPage = 16, byChain = TRUE,
myTheme = theme(text = element_text(size = 12),
axis.text = element_text(size = 10)))
plot_mcmcDensity <- mcmcDensity(popParArray, nParPerPage = 16,
myTheme = theme(text = element_text(size = 12),
axis.text = element_text(size = 10)))
plot_mcmcHistory. [[1]]
Pediatric atorvastatin PK example: MCMC trace plots.
.
. [[2]]
Pediatric atorvastatin PK example: MCMC trace plots.
plot_mcmcDensityByChain. [[1]]
Pediatric atorvastatin PK example: Univariate marginal densities by chain.
plot_mcmcDensity. [[1]]
Pediatric atorvastatin PK example: Univariate marginal densities.
mcmc_pairs(popParArray[,,!grepl("^OMEGA", dimnames(popParArray)[[3]])])
Pediatric atorvastatin PK example: Bivariate marginal densities.
mcmc_pairs(popParArray[,,grepl("^OMEGA", dimnames(popParArray)[[3]])])
Pediatric atorvastatin PK example: Bivariate marginal densities.
caption <- paste("Pediatric atorvastatin PK example:",
c(rep("MCMC trace plots.", length(plot_mcmcHistory)),
rep("Univariate marginal densities by chain.",
length(plot_mcmcDensityByChain)),
rep("Univariate marginal densities.",
length(plot_mcmcDensityByChain)),
rep("Bivariate marginal densities.", 2)))
ptable <- monitor(popParArray,
warmup = 0, print = FALSE) %>%
as.matrix %>%
formatC(3) %>%
as.data.frame
write.csv(ptable, file = file.path(tabDir, paste(modelName,
"ParameterTable.csv", sep = "")))
ptable %>%
rename(SEmean = se_mean, SD = sd, pct2.5 = "2.5%", pct25 = "25%", median = "50%",
pct75 = "75%", pct97.5 = "97.5%", Neff = "n_eff") %>%
mutate(parameter = rownames(.), "95% CI" = paste("(", pct2.5, ", ", pct97.5, ")",
sep = "")) %>%
select(parameter, mean, SD, median, "95% CI", Neff, Rhat) %>%
kable(caption = "Summary of model parameter estimates.") %>%
kable_styling(bootstrap_options = "striped", full_width = F)| parameter | mean | SD | median | 95% CI | Neff | Rhat |
|---|---|---|---|---|---|---|
| OMEGA11 | 0.416 | 0.171 | 0.419 | (0.195, 0.769) | 3 | 1.56 |
| OMEGA21 | 0.0898 | 0.128 | 0.125 | (-0.213, 0.305) | 37 | 1.11 |
| OMEGA22 | 294 | 508 | 1.04 | (0.48, 1.17e+03) | 2 | 1.53 |
| CLHat | 590 | 151 | 618 | ( 373, 842) | 3 | 1.57 |
| V2Hat | 2.44e+03 | 2.36e+03 | 1.17e+03 | ( 667, 6.52e+03) | 2 | 1.53 |
| QHat | 394 | 213 | 293 | ( 188, 755) | 2 | 1.52 |
| V3Hat | 1.6e+03 | 200 | 1.6e+03 | (1.37e+03, 2.01e+03) | 4 | 1.33 |
| KAHat | 0.27 | 0.0396 | 0.281 | (0.209, 0.329) | 2 | 1.53 |
| F1Chewable | 0.896 | 0.179 | 0.87 | (0.675, 1.21) | 374 | 1.55 |
| sigma1 | 0.258 | 0.0494 | 0.235 | (0.208, 0.341) | 2 | 1.53 |
Construct design data frame
## For now just use the NONMEM data set
design <- nmData %>%
filter(C == ".")mrgsolve model code for posterior predictions
library(mrgsolve)
indCode =
'
$PARAM
THETA1 = 5.65
THETA2 = 5.60
THETA3 = 4.29
THETA4 = -1.33
THETA5 = 7.07
THETA6 = 1
THETA7 = -3.2189
WT = 70
FORM = 1
ETA1 = 0
ETA2 = 0
$MAIN
double VWT = log(WT/70);
double CL = exp(THETA1 + (VWT * 0.75) + ETA1 + ETA(1));
double V2 = exp(THETA2 + VWT + ETA2 + ETA(2));
double Q = exp(THETA3 + (VWT * 0.75));
double V3 = exp(THETA5 + VWT);
//double T1 = CL/V2;
//double T23 = Q/V2;
//double T32 = Q/V3;
//double L2 = ((T1+T23+T32)-sqrt(pow(T1+T23+T32, 2)-4*T1*T32))/2;
double KA = exp(THETA4);
F_GUT = 1.0; // TABLET AS REFERENCE
if(FORM == 2) F_GUT = exp(THETA6); // CHEWABLE F1
$PKMODEL cmt = "GUT CENT PERIPH", depot = TRUE
$SIGMA 1
$OMEGA @block
0 0 0
$TABLE
double cPredHat = 1000 * (CENT / V2);
double cPred = cPredHat * (1 + exp(THETA7) * EPS(1));
$CAPTURE cPred
'
## Compile model
atvMod <- mcode("ATVmodel", indCode)
if(FALSE){
## Test runs
out <- atvMod %>%
data_set(design) %>%
carry_out(DOSE) %>%
mrgsim
plot(out, cPred ~ TIME | factor(DOSE), scales = "same")
out <- atvMod %>%
data_set(design) %>%
param(popParameters2[1,] %>%
rename(CL = CLHat,
V2 = V2Hat,
Q = QHat,
V3 = V3Hat,
KA = KAHat,
F1 = F1Chewable)) %>%
carry_out(DOSE) %>%
mrgsim
plot(out, cPred ~ TIME | factor(DOSE), scales = "same")
}Individual predictions
samples <- unique(indParameters$sample)
samples <- 1:8
indSim1 <- mclapply(samples,
function(thisSample){
atvMod %>%
data_set(design) %>%
idata_set(indParameters %>%
filter(sample == thisSample)) %>%
carry_out(sample, ID, TIME, DV, EVID, DOSE, FORM) %>%
mrgsim %>%
as.tbl
}) %>%
bind_rowsPopulation predictions
samples <- unique(popParameters$sample)
popSim1 <- mclapply(samples,
function(thisSample){
atvMod %>%
data_set(design) %>%
param(popParameters %>%
filter(sample == thisSample)) %>%
omat(matrix(with(popParameters %>%
filter(sample == thisSample),
c(OMEGA11, OMEGA21, OMEGA21, OMEGA22)), ncol = 2)) %>%
carry_out(sample, ID, TIME, DV, EVID, DOSE, FORM) %>%
mrgsim %>%
as.tbl
}) %>%
bind_rowsPosterior predictive distributions
## Combined plots
indPred <- indSim1 %>%
group_by(ID, TIME, EVID, DOSE, FORM) %>%
summarize(lbInd = quantile(cPred, probs = 0.05, na.rm = TRUE),
medianInd = quantile(cPred, probs = 0.5, na.rm = TRUE),
ubInd = quantile(cPred, probs = 0.95, na.rm = TRUE))
popPred <- popSim1 %>%
group_by(ID, TIME, EVID, DOSE, FORM) %>%
summarize(lbPop = quantile(cPred, probs = 0.05, na.rm = TRUE),
medianPop = quantile(cPred, probs = 0.5, na.rm = TRUE),
ubPop = quantile(cPred, probs = 0.95, na.rm = TRUE))
predAll <- indPred %>%
left_join(popPred) %>%
left_join(design) %>%
mutate(DV = suppressWarnings(as.numeric(DV)))
forms <- sort(unique(predAll$FORM))
plot_PPCPK <- lapply(forms,
function(thisForm){
p1 <- ggplot(predAll %>% filter(FORM == thisForm),
aes(x = TIME, y = DV))
p1 <- p1 +
geom_line(aes(x = TIME, y = medianPop,
color = "population")) +
geom_ribbon(aes(ymin = lbPop, ymax = ubPop,
fill = "population"),
alpha = 0.25) +
geom_line(aes(x = TIME, y = medianInd,
color = "individual")) +
geom_ribbon(aes(ymin = lbInd, ymax = ubInd,
fill = "individual"),
alpha = 0.25) +
scale_color_brewer(name ="prediction",
breaks=c("individual", "population"),
palette = "Set1") +
scale_fill_brewer(name ="prediction",
breaks=c("individual", "population"),
palette = "Set1")
p1 + geom_point() +
labs(title = paste("formulation =", c("tablet", "chewable")[thisForm]),
x = "time (h)",
y = "atorvastatin plasma concentration (nM)") +
theme(text = element_text(size = 12),
axis.text = element_text(size = 12),
## legend.position = c(0.8, 0.25),
strip.text = element_text(size = 8)) +
facet_wrap(~ ID, scales = "free_y")
})
plot_PPCPK[[1]]
. Warning: Removed 48 rows containing missing values (geom_path).. Warning: Removed 367 rows containing missing values (geom_point).
[[2]]
. Warning: Removed 227 rows containing missing values (geom_point).
sessionInfo(). R version 3.5.1 (2018-07-02)
. Platform: x86_64-pc-linux-gnu (64-bit)
. Running under: Ubuntu 14.04.5 LTS
.
. Matrix products: default
. BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
. LAPACK: /usr/lib/lapack/liblapack.so.3.0
.
. locale:
. [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
. [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
. [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
. [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
. [9] LC_ADDRESS=C LC_TELEPHONE=C
. [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
.
. attached base packages:
. [1] parallel grid stats graphics grDevices utils datasets
. [8] methods base
.
. other attached packages:
. [1] mrgsolve_0.9.2 PKPDmisc_2.1.1 qapply_1.39
. [4] rlecuyer_0.3-4 kableExtra_1.1.0 knitr_1.25
. [7] gridExtra_2.3 forcats_0.4.0 stringr_1.4.0
. [10] dplyr_0.8.3 purrr_0.3.2 readr_1.3.1
. [13] tidyr_1.0.0 tibble_2.1.3 tidyverse_1.2.1
. [16] bayesplot_1.7.0 rstan_2.19.2 ggplot2_3.2.1
. [19] StanHeaders_2.19.0 metrumrg_5.57 MASS_7.3-51.1
. [22] XML_3.98-1.20 lattice_0.20-38 reshape_0.8.8
.
. loaded via a namespace (and not attached):
. [1] httr_1.4.1 jsonlite_1.6
. [3] viridisLite_0.3.0 modelr_0.1.5
. [5] assertthat_0.2.1 highr_0.8
. [7] stats4_3.5.1 cellranger_1.1.0
. [9] yaml_2.2.0 pillar_1.4.2
. [11] backports_1.1.5 glue_1.3.1
. [13] digest_0.6.21 RColorBrewer_1.1-2
. [15] rvest_0.3.4 colorspace_1.4-1
. [17] htmltools_0.4.0 plyr_1.8.4
. [19] pkgconfig_2.0.3 broom_0.5.2
. [21] haven_2.1.1 scales_1.0.0
. [23] webshot_0.5.1 processx_3.4.1
. [25] generics_0.0.2 withr_2.1.2
. [27] lazyeval_0.2.2 cli_1.1.0
. [29] magrittr_1.5 crayon_1.3.4
. [31] readxl_1.3.1 evaluate_0.14
. [33] ps_1.3.0 nlme_3.1-137
. [35] RcppArmadillo_0.9.800.1.0 xml2_1.2.2
. [37] pkgbuild_1.0.6 data.table_1.12.4
. [39] tools_3.5.1 loo_2.1.0
. [41] prettyunits_1.0.2 hms_0.5.1
. [43] lifecycle_0.1.0 matrixStats_0.55.0
. [45] munsell_0.5.0 callr_3.3.2
. [47] compiler_3.5.1 rlang_0.4.0
. [49] ggridges_0.5.1 rstudioapi_0.10
. [51] labeling_0.3 rmarkdown_1.16
. [53] gtable_0.3.0 inline_0.3.15
. [55] reshape2_1.4.3 R6_2.4.0
. [57] lubridate_1.7.4 zeallot_0.1.0
. [59] stringi_1.4.3 Rcpp_1.0.2
. [61] vctrs_0.2.0 tidyselect_0.2.5
. [63] xfun_0.10